library(Seurat)
library(tidyseurat)
library(tidyverse)
library(ggpubr)
library(patchwork)
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

test_div_percent <- function(xina, yina, xinb, yinb){
  seq1 <- c(rep(1,xina),rep(0,yina))
  seq2 <- c(rep(1,xinb),rep(0,yinb))
  conf1 <- t.test(seq1) |>
    broom::tidy() |>
    select(estimate, conf.low, conf.high)
  conf2 <- t.test(seq2) |>
    broom::tidy() |>
    select(estimate, conf.low, conf.high)
  
  bind_rows(conf1, conf2, .id = 'id') |>
    mutate(pval = t.test(seq1, seq2)$p.value)
}

test_for3 <- function(df){
  res <- df |>
    pivot_wider(names_from = genotype, values_from = c(n,other)) |>
    reframe(ab = test_div_percent(n_IT,other_IT,n_II,other_II),
            ac = test_div_percent(n_TT,other_TT,n_II,other_II),
            bc = test_div_percent(n_IT,other_IT,n_TT,other_TT)) |>
    map(slice_max, conf.high) |>
    list_rbind() |>
    mutate(group1 = c('II','II','IT'),
           group2 = c('IT','TT','TT'),
           y.position = max(conf.high) * 1.1,
           pval = case_when(pval < .001 ~ '***',
                            pval < .01 ~ '**',
                            pval < .05 ~ '*', .default = 'ns'),
           .keep = 'none')
}

i2t_colors <- c('forestgreen','royalblue','red2')

meta.c4 <- read_csv('CRC-I/results/crc_merge4_meta.csv')
meta.myl4 <- read_csv('CRC-I/results/crc_merge4_myel.meta.csv')
meta.bc4 <- read_csv('CRC-I/results/crc_merge4_bc.meta.csv')
meta.nkt4 <- read_csv('CRC-I/results/crc_merge4_nkt.meta.csv')

fcgr2b.posi <- read_csv('CRC-I/results/merge4.fcgr2b.pos.cell.csv')

geno.c4 <- meta.c4 |>
  select(.cell, genotype)

# tidy and merge meta ----------
plot_frac <- function(df, cell.type){
  df |>
    calc_frac_conf_on_grouped_count(.data$genotype, {{ cell.type }}) |>
    ggplot(aes(.data$genotype, .data$fraction, fill = .data$genotype,
               ymax = .data$conf.high, ymin = .data$conf.low)) +
    geom_col() +
    geom_errorbar(width = .3) +
    facet_wrap(vars({{ cell.type }}), scales = 'free_y') +
    theme_pubr() +
    scale_fill_manual(values = i2t_colors)
}

## major types --------
meta.c4 |>
  filter(.cell %in% fcgr2b.posi$.cell) |>
  calc_frac_conf_on_grouped_count(genotype, manual.main) |>
  ggplot(aes(genotype, fraction,
             ymax = conf.high, ymin = conf.low, fill = genotype)) +
  geom_col() +
  geom_errorbar(width = .3) +
  facet_wrap(~manual.main, scales = 'free_y') +
  theme_pubr() +
  scale_fill_manual(values = i2t_colors) +
  labs(title = 'Major immune cells (FCGR2B+)')

## myeloid ------------
meta.myl4 |>
  left_join(geno.c4) |>
  filter(.cell %in% fcgr2b.posi$.cell) |>
  plot_frac(namely.fine) +
  ggtitle('Major myeloid cells (FCGR2B+)')

meta.myl4 |>
  left_join(geno.c4) |>
  filter(.cell %in% fcgr2b.posi$.cell) |>
  plot_frac(zhang.fine) +
  ggtitle('Subtypes of myeloid cells (FCGR2B+)')

## T cell ---------
meta.nkt4 |>
  left_join(geno.c4) |>
  filter(.cell %in% fcgr2b.posi$.cell) |>
  plot_frac(namely.fine) +
  ggtitle('Major T/NK cells (FCGR2B+)')

meta.nkt4 |>
  filter(str_detect(zhang.fine, 'CD4'),
         .cell %in% fcgr2b.posi$.cell) |>
  left_join(geno.c4) |>
  plot_frac(zhang.fine) +
  ggtitle('Subtypes of CD4+ T cells (FCGR2B+)')

meta.nkt4 |>
  filter(!str_detect(zhang.fine, 'CD4'),
         .cell %in% fcgr2b.posi$.cell) |>
  left_join(geno.c4) |>
  plot_frac(zhang.fine) +
  ggtitle('Subtypes of CD8+ T/NK cells (FCGR2B+)')

## B cell ---------
meta.bc4 |>
  filter(.cell %in% fcgr2b.posi$.cell) |>
  left_join(geno.c4) |>
  plot_frac(namely.fine) +
  ggtitle('Major B cells (FCGR2B+)')

meta.bc4 |>
  left_join(geno.c4) |>
  filter(.cell %in% fcgr2b.posi$.cell) |>
  plot_frac(zhang.fine) +
  ggtitle('Subtypes of B cells (FCGR2B+)')

## CCL22+ DC3 ---------
sg_sobj_dc <- sobj_sg |>
  filter(str_detect(zhang2020_fine, 'DC'))

dc_fine <- sg_sobj_dc |>
  distinct(seurat_clusters) |>
  arrange(seurat_clusters) |>
  mutate(dc_cluster = seq_along(seurat_clusters),
         dc_cluster = str_glue('DC cluster {dc_cluster}'))

sg_sobj_dc <- sg_sobj_dc |>
  left_join(dc_fine)

sg_sobj_dc |>
  VlnPlot('CCL22', group.by = 'dc_cluster')

sg_sobj_dc |>
  count(dc_cluster, zhang2020_fine)

dc3_marker <- sg_sobj_dc |>
  FindMarkers(ident.1 = 'DC cluster 5', group.by = 'dc_cluster') |>
  as_tibble(rownames = 'gene')

dc3_marker |>
  EnhancedVolcano::EnhancedVolcano(lab = dc3_marker$gene, x = 'avg_log2FC', y = 'p_val_adj') +
  labs(title = 'CCL22+ DC3 differential expressed genes vs other DCs',
       subtitle = NULL)

ggplot(aes(avg_log2FC,-log10(p_val_adj))) +
  geom_point()

updc3 <- dc3_marker |>
  dplyr::filter(avg_log2FC > 0, p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           minGSSize = 3,
           readable = TRUE)

updc3@result |> 
  as_tibble() |>
  head(n = 10) |>
  mutate(Description = fct_reorder(Description, Count)) |>
  ggplot(aes(Description, Count, fill = p.adjust)) +
  geom_col() +
  coord_flip() +
  theme_pubr() +
  labs_pubr() +
  scale_fill_gradient(low = 'red', high = 'black') +
  ggtitle('Upregulated pathways in CCL22+ DC3')

sobj_dc <- sobj |>
  filter(str_detect(Sub_Cluster, 'DC'))

sobj_dc <- sobj_dc |>
  mutate(orig.ident = as.character(orig.ident)) |>
  FindVariableFeatures() |>
  ScaleData() |>
  RunPCA() |>
  harmony::RunHarmony('orig.ident') |>
  RunUMAP(reduction = "harmony", dims = 1:20) |>
  RunTSNE(reduction = "harmony", dims = 1:20) |>
  FindNeighbors(reduction = "harmony", dims = 1:25) |>
  FindClusters(resolution = 2)

sobj_dc |> VlnPlot('CCL22')

sobj_dc |> count(RNA_snn_res.2,Sub_Cluster)

ccl22_id10x <- sobj_dc |> filter(RNA_snn_res.2 == 1) |> pull(.cell)

sg_sobj_dc |> VlnPlot('CCL22',group.by = 'RNA_snn_res.2')

sg_sobj_dc |> count(RNA_snn_res.2)

ccl22_idsg <- sg_sobj_dc |>
  filter(RNA_snn_res.2 == 18) |> pull(.cell)

meta_all4 <- meta_all4 |>
  mutate(type_fine = ifelse(ID %in% c(ccl22_id10x,ccl22_idsg),
                            'DC3_CCL22', type_fine))

meta_all4 |>
  write_csv('CRC-I/results/crc-cell-type-all4.csv')

# meta file checkpoint ------------
meta_all4 <- read_csv('CRC-I/results/crc-cell-type-all4.csv')

meta_pub4 <- read_csv('CRC-I/results/crc_merge4_meta.csv')

meta_pub4 |> ggplot(aes(genotype, fill = genotype)) +
  geom_bar() +
  theme_pubr() +
  labs_pubr() +
  scale_fill_manual(values = i2t_colors)

meta_pub4 |> write_csv('CRC-I/results/crc_merge4_meta.csv')

meta_all4 |> count(type_main)

meta_imm4 <- meta_all4 |>
  mutate(type_main = case_when(type_main %in% c('B cell', 'Myeloid cell') ~ type_main,
                               str_detect(type_main, 'CD4') ~ 'CD4 T cell',
                               str_detect(type_main, 'CD8') ~ 'CD8 T cell',
                               .default = 'NA')) |>
  filter(type_main != 'NA' & !is.na(type_fine))

meta_imm4 |> count(type_main)

meta_imm4 |>
  filter(type_main == 'CD8 T cell') |>
  count(type_fine, sort = T)

meta_imm3 <- meta_imm4 |>
  filter(Platform != 'smart2018')

patient2genotype <- meta_imm3 |>
  count(patient, genotype) |>
  select(-n)

# major cell type by patient ---------
major_type_patient <-meta_imm3 |>
  group_by(patient) |>
  count(type_main) |>
  mutate(sum = sum(n), frac = n/sum) |>
  left_join(patient2genotype)

major_type_patient |>
  ggplot(aes(genotype, frac)) +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~type_main)
# it won't work: some patient have too few cells to have signficance

# only 7 patient have 1000+ cells: 4 II + 1 IT + 2 TT
sig_pat <- meta_imm4 |>
  count(patient, Platform, genotype) |>
  filter(n > 500) |>
  pull(patient)

major_type_patient |>
  filter(patient %in% sig_pat) |>
  ggplot(aes(genotype, frac, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(method = 't.test',
                     comparisons = list(c('II','TT'),c('IT','TT'))) +
  facet_wrap(~type_main, scales = 'free') +
  theme_pubr() +
  scale_y_continuous(expand = expansion(mult = c(.1,.2))) +
  scale_color_manual(values = c('green3','blue','red')) +
  labs(y = 'fraction in TME immune cells')

# cell frac in conf interval ------
calc_conf <- function(df){
  df |>
    mutate(other = sum(n) - n,
           fraction = n / sum(n)) |>
    rowwise() |>
    mutate(conf.low = t.test(c(rep(1,n),rep(0,other)))$conf.int[1],
           conf.high = t.test(c(rep(1,n),rep(0,other)))$conf.int[2])
}

meta_conf <- meta_imm3 |>
  group_by(genotype) |>
  count(type_main) |>
  mutate(other = sum(n) - n,
         fraction = n / sum(n)) |>
  rowwise() |>
  mutate(conf.low = t.test(c(rep(1,n),rep(0,other)))$conf.int[1],
         conf.high = t.test(c(rep(1,n),rep(0,other)))$conf.int[2])

pub4_conf <- meta_pub4 |>
  group_by(genotype) |>
  dplyr::count(manual_main) |>
  calc_conf()

calc_pval <- function(df){
  df |>
    mutate(other = sum(n) - n) |>
    nest(data = c(genotype, n, other)) |>
    rowwise() |>
    mutate(data = list(test_for3(data))) |>
    unnest(data)
}

pub4_pval <- meta_pub4 |>
  group_by(genotype)|>
  dplyr::count(manual_main) |>
  calc_pval()

conf_pval <- meta_imm3 |>
  group_by(genotype) |>
  count(type_main) |>
  mutate(other = sum(n) - n) |>
  nest(data = c(genotype,n)) |>
  rowwise() |>
  mutate(data = list(test_for3(data))) |>
  unnest(data)

pub4_conf |>
  filter((manual_main %in% c('Mast','Myeloid','NK'))) |>
  ggplot(aes(genotype, fraction,
             ymin = conf.low, ymax = conf.high, color = genotype)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .3) +
  facet_wrap(~manual_main, scales = 'free')+
  theme_pubr() +
  labs_pubr() +
  stat_pvalue_manual(data = filter(pub4_pval, manual_main %in% c('Mast','Myeloid','NK')),
                     label = 'pval',hide.ns = 'pval',
                     step.increase = .2,
                     step.group.by = 'manual_main',
                     inherit.aes = FALSE) +
  scale_y_continuous(expand = expansion(mult = c(0,.2))) +
  scale_color_manual(values = c('green3','blue','red')) +
  labs(y = 'fraction in TME immune cells')

# subtype in B cell ---------
b2_conf <- meta_imm3 |>
  filter(type_main == 'B cell') |>
  group_by(genotype) |>
  count(type_fine) |>
  mutate(other = sum(n) - n,
         fraction = n / sum(n)) |>
  rowwise() |>
  mutate(conf.low = t.test(c(rep(1,n),rep(0,other)))$conf.int[1],
         conf.high = t.test(c(rep(1,n),rep(0,other)))$conf.int[2]) |>
  filter(str_detect(type_fine, 'IgA|IgG|MS4A1'))


b2_pval <- meta_imm3 |>
  filter(type_main == 'B cell') |>
  group_by(genotype) |>
  count(type_fine) |>
  mutate(other = sum(n) - n) |>
  filter(str_detect(type_fine, 'IgA|IgG|MS4A1')) |>
  nest(data = -type_fine) |>
  rowwise() |>
  mutate(data = list(test_for3(data))) |>
  unnest(data)

b2_conf |>
  ggplot(aes(genotype, fraction,
             ymin = conf.low, ymax = conf.high, color = genotype)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .3) +
  facet_wrap(~type_fine, scales = 'free')+
  theme_pubr() +
  labs_pubr() +
  stat_pvalue_manual(data = b2_pval,
                     label = 'pval',hide.ns = 'pval',
                     step.increase = .2,
                     step.group.by = 'type_fine',
                     inherit.aes = FALSE) +
  scale_y_continuous(expand = expansion(mult = c(0,.2))) +
  scale_color_manual(values = c('green3','blue','red')) +
  labs(y = 'Fraction in all TME B cells')

# subtype in myeloid cells -----
## by namely subtypes ------
meta_myl <- meta_imm3 |>
  filter(str_detect(type_main, 'yeloid')) |>
  mutate(type_namely = str_extract(type_fine, 'Macro|Mast|Mono|TAM|cDC|pDC'),
         type_namely = case_match(type_namely, 'TAM' ~ 'Macro',
                                  .default = type_namely))

meta_myl |>
  ggplot(aes(genotype, fill = genotype)) +
  geom_bar() +
  scale_fill_manual(values = i2t_colors) +
  facet_wrap(~ type_namely, scales = 'free_y') +
  theme_pubr() +
  labs(title = 'Count of myeloid in CRC TME')

meta_myl |> write_csv('CRC-I/results/crc_merge3_myl.namely.csv')

## by snn clusters ---------
pub4_myl_conf <- meta_pub4 |>
  filter(manual_main == 'Myeloid') |>
  group_by(genotype) |>
  dplyr::count(seurat_clusters) |>
  calc_conf()

pub4_myl_pval <- meta_pub4 |>
  filter(manual_main == 'Myeloid') |>
  group_by(genotype) |>
  dplyr::count(seurat_clusters) |>
  calc_pval()

pub4_myl_conf |>
  ggplot(aes(genotype, fraction,
             ymin = conf.low, ymax = conf.high, color = genotype)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .3) +
  facet_wrap(~seurat_clusters, scales = 'free', ncol = 4)+
  theme_pubr() +
  labs_pubr() +
  stat_pvalue_manual(data = pub4_myl_pval,
                     label = 'pval',hide.ns = 'pval',
                     step.increase = .2,
                     step.group.by = 'seurat_clusters',
                     inherit.aes = FALSE) +
  scale_y_continuous(expand = expansion(mult = c(0,.2))) +
  scale_color_manual(values = i2t_colors) +
  labs(y = 'Fraction in all TME myeloid cells')

sig_submyl <- meta_imm3 |>
  filter(type_main == 'Myeloid cell') |>
  count(type_fine, genotype) |>
  count(type_fine) |>
  filter(n == 3) |>
  pull(type_fine)

myl_pval <- meta_imm3 |>
  filter(type_main == 'Myeloid cell') |>
  count(type_fine, genotype) |>
  group_by(genotype) |>
  mutate(other = sum(n) - n) |>
  filter(type_fine %in% sig_submyl) |>
  nest(data = -type_fine) |>
  rowwise() |>
  mutate(data = list(test_for3(data))) |>
  unnest(data)

myl_conf |>
  filter(type_fine %in% sig_submyl) |>
  ggplot(aes(genotype, fraction,
             ymin = conf.low, ymax = conf.high, color = genotype)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .3) +
  facet_wrap(~type_fine, scales = 'free', ncol = 4)+
  theme_pubr() +
  labs_pubr() +
  stat_pvalue_manual(data = myl_pval,
                     label = 'pval',hide.ns = 'pval',
                     step.increase = .2,
                     step.group.by = 'type_fine',
                     inherit.aes = FALSE) +
  scale_y_continuous(expand = expansion(mult = c(0,.2))) +
  scale_color_manual(values = c('green3','blue','red')) +
  labs(y = 'Fraction in all TME myeloid cells')

# subtype in CD4 cells -----
cd4_imm3 <- meta_imm4 |>
  filter(type_main == 'CD4 T cell')

sig_cd4 <- cd4_imm3 |>
  count(type_fine, genotype) |>
  count(type_fine) |>
  filter(n == 3) |>
  pull(type_fine)

cd4_conf <- cd4_imm3 |>
  group_by(genotype) |>
  count(type_fine) |>
  mutate(other = sum(n) - n,
         fraction = n / sum(n)) |>
  rowwise() |>
  mutate(conf.low = t.test(c(rep(1,n),rep(0,other)))$conf.int[1],
         conf.high = t.test(c(rep(1,n),rep(0,other)))$conf.int[2])

cd4_pval <- cd4_imm3 |>
  count(type_fine, genotype) |>
  group_by(genotype) |>
  mutate(other = sum(n) - n) |>
  filter(type_fine %in% sig_cd4) |>
  nest(data = -type_fine) |>
  rowwise() |>
  mutate(data = list(test_for3(data))) |>
  unnest(data)

cd4_conf |>
  filter(type_fine %in% sig_cd4) |>
  ggplot(aes(genotype, fraction,
             ymin = conf.low, ymax = conf.high, color = genotype)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .3) +
  facet_wrap(~type_fine, scales = 'free', ncol = 4)+
  theme_pubr() +
  labs_pubr() +
  stat_pvalue_manual(data = cd4_pval,
                     label = 'pval',hide.ns = 'pval',
                     step.increase = .2,
                     step.group.by = 'type_fine',
                     inherit.aes = FALSE) +
  scale_y_continuous(expand = expansion(mult = c(0,.2))) +
  coord_cartesian(ylim = c(0,NA)) +
  scale_color_manual(values = c('green3','blue','red')) +
  labs(y = 'Fraction in TME CD4+ cells')

# subtype in CD8 ---------
cd8_imm3 <- meta_imm4 |>
  filter(type_main == 'CD8 T cell')

sig_cd8 <- cd8_imm3 |>
  count(type_fine, genotype) |>
  count(type_fine) |>
  filter(n == 3) |>
  pull(type_fine)

cd8_conf <- cd8_imm3 |>
  group_by(genotype) |>
  count(type_fine) |>
  mutate(other = sum(n) - n,
         fraction = n / sum(n)) |>
  rowwise() |>
  mutate(conf.low = t.test(c(rep(1,n),rep(0,other)))$conf.int[1],
         conf.high = t.test(c(rep(1,n),rep(0,other)))$conf.int[2])

cd8_pval <- cd8_imm3 |>
  count(type_fine, genotype) |>
  group_by(genotype) |>
  mutate(other = sum(n) - n) |>
  filter(type_fine %in% sig_cd8) |>
  nest(data = -type_fine) |>
  rowwise() |>
  mutate(data = list(test_for3(data))) |>
  unnest(data)

cd8_conf |>
  filter(type_fine %in% sig_cd8) |>
  ggplot(aes(genotype, fraction,
             ymin = conf.low, ymax = conf.high, color = genotype)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .3) +
  facet_wrap(~type_fine, scales = 'free', ncol = 3)+
  theme_pubr() +
  labs_pubr() +
  stat_pvalue_manual(data = cd8_pval,
                     label = 'pval',hide.ns = 'pval',
                     step.increase = .2,
                     step.group.by = 'type_fine',
                     inherit.aes = FALSE) +
  scale_y_continuous(expand = expansion(mult = c(0,.2))) +
  scale_color_manual(values = c('green3','blue','red')) +
  coord_cartesian(ylim = c(0,NA)) +
  labs(y = 'Fraction in TME CD8+ cells')

# CD56dim CD16hi DNAJB1+ TaNK ---------
tank_sg <- sobj |>
  get_abundance_sc_wide(c('NCAM1','KLRF1','DNAJB1','HSPA1A','CD3D','CD3G','CD3E')) |>
  filter(NCAM1 & KLRF1 & DNAJB1 & CD3E + CD3G + CD3D == 0) |>
  pull(.cell)

tank_10x <- sobj_10x |>
  get_abundance_sc_wide(c('NCAM1','KLRF1','DNAJB1','HSPA1A','CD3D','CD3G','CD3E')) |>
  filter(NCAM1 & KLRF1 & DNAJB1 & CD3E + CD3G + CD3D == 0) |>
  pull(.cell)

tank_sm2020 <- sobj_sm2020 |>
  get_abundance_sc_wide(c('NCAM1','KLRF1','DNAJB1','HSPA1A','CD3D','CD3G','CD3E')) |>
  filter(NCAM1 & KLRF1 & DNAJB1 & CD3E + CD3G + CD3D == 0 & str_detect(.cell, 'T')) |>
  pull(.cell)

sobj_10x |> count(orig.ident, genotype)
sobj_sm2020 |> count(orig.ident, genotype)

sobj_sg |> filter(str_detect(zhang2020_main, 'ILC')) |> count(genotype)
sobj_10x |> filter(str_detect(Sub_Cluster, 'NK')) |> count(genotype)
sobj_sm2020 |> filter(str_detect(Sub_Cluster, 'NK') & Tissue == 'T') |> count(genotype)

tank_meta <- tibble(genotype = c('II','IT','TT'), n = c(9,10,1), other = c(208-9,166-10,129-1))

pval_tank <- tank_meta |> test_for3()

conf_tank <- tank_meta |>
  rowwise() |>
  mutate(conf.low = t.test(c(rep(1,n),rep(0,other)))$conf.int[1],
         conf.high = t.test(c(rep(1,n),rep(0,other)))$conf.int[2])

conf_tank |>
  mutate(fraction = n / (n + other)) |>
  ggplot(aes(genotype, fraction,
             ymin = conf.low, ymax = conf.high, color = genotype)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .3) +
  theme_pubr() +
  labs_pubr() +
  stat_pvalue_manual(data = pval_tank,
                     label = 'pval',hide.ns = 'pval',
                     step.increase = .2,
                     inherit.aes = FALSE) +
  scale_y_continuous(expand = expansion(mult = c(0,.2))) +
  scale_color_manual(values = c('green3','blue','red')) +
  coord_cartesian(ylim = c(0,NA)) +
  labs(y = 'Fraction in TME NK cells')

g1 <- last_plot()

conf_tank |>
  ggplot(aes(genotype, n, fill = genotype)) +
  geom_col() +
  geom_text(aes(label = n), color = 'white', vjust = 2, size = 6) +
  theme_pubr() +
  labs_pubr() +
  scale_y_continuous(breaks = c(0,5,10)) +
  scale_fill_manual(values = c('green3','blue','red')) +
  labs(y = 'Number of CD56dim CD16hi DNAJB1+ TaNK cells')

g2 <- last_plot()

g1 + g2 + plot_annotation(title = 'CD56dim CD16hi DNAJB1+ TaNK cells in CRC')

# join with ecotype states --------
recovery_path <- c('../ecotyper/Recovery_seek3v2',
                   '../ecotyper/Recovery_zhang2020_10x/',
                   '../ecotyper/Recovery_zhang2018_smart/')

eco_meta <- list.files(recovery_path,
                       pattern = 'state_assignment.txt',
                       full.names = TRUE,
                       recursive = TRUE) |>
  read_delim() |>
  select(-InitialState)

meta_imm3 <- meta_imm3 |>
  left_join(eco_meta)

meta_imm3 |>
  filter(type_main == 'B cell') |>
  ggplot(aes(fill =State, x = type_fine)) +
  geom_bar(position = 'fill')

write_csv(meta_imm3, 'CRC-I/results/sg+zhang_imm3_meta.csv')

## cell-states in cd8 -----------
cd8_imm3 <- meta_imm3 |>
  filter(type_main == 'CD8 T cell')

cd8_conf <- cd8_imm3 |>
  group_by(genotype) |>
  count(State) |>
  mutate(other = sum(n) - n,
         fraction = n / sum(n)) |>
  rowwise() |>
  filter(!is.na(State)) |>
  mutate(conf.low = t.test(c(rep(1,n),rep(0,other)))$conf.int[1],
         conf.high = t.test(c(rep(1,n),rep(0,other)))$conf.int[2])

cd8_pval <- cd8_imm3 |>
  count(State, genotype) |>
  group_by(genotype) |>
  mutate(other = sum(n) - n) |>
  nest(data = -State) |>
  rowwise() |>
  filter(!is.na(State)) |>
  mutate(data = list(test_for3(data))) |>
  unnest(data)

cd8_conf |>
  ggplot(aes(genotype, fraction,
             ymin = conf.low, ymax = conf.high, color = genotype)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .3) +
  facet_wrap(~State, scales = 'free', ncol = 4)+
  theme_pubr() +
  labs_pubr() +
  stat_pvalue_manual(data = cd8_pval,
                     label = 'pval',hide.ns = 'pval',
                     step.increase = .2,
                     step.group.by = 'State',
                     inherit.aes = FALSE) +
  scale_y_continuous(expand = expansion(mult = c(0,.2))) +
  scale_color_manual(values = c('green3','blue','red')) +
  labs(y = 'Fraction in TME CD8+ cells')

# integrated fraction --------
merge4.meta <- read_csv('CRC-I/results/crc_merge4_meta.csv')

merge4.meta <-
read_csv('CRC-I/results/crc_merge4_mast_meta.csv') |>
  select(.cell, mast.subtype) |>
  right_join(merge4.meta)

## TODO: correct all genotypes in li2021 and guo2021 =========
merge4.meta %<>%
  mutate(genotype = case_match(orig.ident,
                               c('COL07','COL16','COL17') ~ 'IT',
                               c('COL12','COL15','COL18','right1','right2','right3') ~ 'II',
                               .default = genotype))

merge4.meta |>
  calc_frac_conf_on_grouped_count(genotype, hpca_main_c4) |>
  filter(hpca_main_c4 == 'Mast cell') |>
  ggplot(aes(genotype, fraction, fill = genotype)) +
  geom_col() +
  geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width = .3) +
  theme_pubr() +
  scale_fill_manual(values = i2t_colors) +
  labs(title = 'Mast cell in all TME immune cells')

merge4.meta |>
  calc_frac_conf_on_grouped_count(genotype, mast.subtype) |>
  filter(mast.subtype == 'VEGFA+') |>
  ggplot(aes(genotype, fraction, fill = genotype)) +
  geom_col() +
  geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width = .3) +
  theme_pubr() +
  scale_fill_manual(values = i2t_colors) +
  labs(title = 'VEGFA+ Mast cell in all TME immune cells')

merge4.meta |>
  calc_frac_conf_on_grouped_count(genotype, mast.subtype) |>
  filter(str_detect(mast.subtype, '\\+')) |>
  pivot_wider(id_cols = genotype, names_from = mast.subtype, values_from = fraction) |>
  mutate(ratio = `VEGFA+`/`TNF+`) |>
  ggplot(aes(genotype, ratio, fill = genotype)) +
  geom_col() +
  theme_pubr() +
  scale_fill_manual(values = i2t_colors) +
  labs(title = 'Ratio of VEGFA+/TNF+ Mast cell in CRC TME')

## in all myeloid instead of all immune cells
merge4.meta.myl <- merge4.meta |>
  filter(!str_detect(hpca_main_c4, 'B|T|NK'))

merge4.meta.myl |>
  calc_frac_conf_on_grouped_count(genotype, hpca_main_c4) |>
  filter(hpca_main_c4 == 'Mast cell') |>
  ggplot(aes(genotype, fraction, fill = genotype)) +
  geom_col() +
  geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width = .3) +
  theme_pubr() +
  scale_fill_manual(values = i2t_colors) +
  labs(title = 'Mast cell in all TME myeloid cells')

merge4.meta.myl |>
  calc_frac_conf_on_grouped_count(genotype, mast.subtype) |>
  filter(mast.subtype == 'VEGFA+') |>
  ggplot(aes(genotype, fraction, fill = genotype)) +
  geom_col() +
  geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width = .3) +
  theme_pubr() +
  scale_fill_manual(values = i2t_colors) +
  labs(title = 'VEGFA+ Mast cell in all TME myeliod cells')
